/*******************************************************************************
  ARTICLE	GAY, GOBBI, GONI (2025) "REVOLUTIONARY TRANSITIONS. INHERITANCE    
            CHANGE AND FERTILITY DECLINE" JOURNAL OF POLITICAL ECONOMY         
                                                                               
  AUTHORS	VICTOR GAY, PAULA GOBBI, MARC GONI                                 
  CONTACT	victor.gay@tse-fr.eu; paula.eugenia.gobbi@ulb.be; marc.goni@uib.no 
  VERSION	1.0 (MAY 2025)                                                     
  SOFTWARE	STATA SE 18                                                        
  LICENCE	MIT                                                                
--------------------------------------------------------------------------------

GENI RESULTS PAPER DO FILE

This file contains the codes to reproduce the figures and tables in the paper using the Geni data.

Instructions: 
-------------
	Gain access to the Geni database from MyHeritage, Ltd and either:
	
	1) place the following files:

		geni_profiles.csv (20.2 GB, created on April 4, 2022, at 11:39:31),
		geni_unions.csv (7.2 GB, created on April 4, 2022, at 11:55:59),
		geni_union details.csv (2.6 GB, created on April 4, 2022, at 11:54:34).
	
	in folder /1_raw_data/1_1_henri/ (see README for more details) and run R-codes 
	named 1-* to 8-* in folder "/2_scripts/2_1_data/01_geni_data_to_sample"; or
	
	2) place the author-provided fr-clean.csv file into the folder 3_outputs/3_1_datasets.
 
	After 1) or 2), run do-file 02_geni-inheritance-and-controls.do and 03_geni-final-data-prep.do	
	
	Open do-files from directory where they are placed; order matters; run whole code.
	
Contents: 
---------
	Program setup

	A. Supplementary tables
		Table A9: Heterogeneous effects by soil conditions for small versus large farms, RD-DD, Geni database.

	E. Robustness of RD-DD results, Geni database
		Table E1: Covariate balance test for spatial regression-discontinuity analysis.
		Table E2: Sensitivity to additional RD-DD specifications.
		Figure E1: Balance RD plots.
		Figure E2: Trends in completed fertility under reformed and not reformed inheritance
		Figure E3: Sensitivity to bandwidth choice for spatial RD-DD.
		Figure E4: Conley adjusted standard errors with different distance cutoffs.

Date last update: May 2025; Ran using STATA 18.5
*/
********************************************************************************


********************
* 0. PROGRAM SETUP *
********************

version 18
clear all
set more off

************************
* PACKAGE DEPENDENCIES *
************************

ssc install rdrobust, replace
ssc install reghdfe, replace
ssc install ftools, replace
ssc install outreg2, replace
ssc install grc1leg2, replace
ssc install coefplot, replace
ssc install acreg, replace
ssc install ranktest, replace
ssc install hdfe, replace

***************
* DIRECTORIES *
***************

global DAT 	= "../../3_outputs/3_1_datasets"
global APPTAB 	= "../../3_outputs/3_3_appendix/3_3_1_appendix_tables"
global APPFIG 	= "../../3_outputs/3_3_appendix/3_3_2_appendix_figures"
global TEMP = "../2_0_tempfiles"

timer on 1

********************************************************************************
*
* APPENDIX A. ADDITIONAL TABLES
*
********************************************************************************


* ==============================================================================
* Table A9: Heterogeneous effects by soil conditions for small versus large farms,
* RD-DD, Geni database.
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)


* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1)-(3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(ruggedness byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (4) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(ruggedness byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(ruggedness byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)


* RD-DD heterogeneity
* -------------------------------
gen TxFxlarge = 0
gen TxFxsmall = 0
label var TxFxlarge "RD-DD estimate (soil conditions for large farms)"
label var TxFxsmall "RD-DD estimate (soil conditions for small farms)"
global treatvars TxFxlarge TxFxsmall

// based on sandy soils
bys insee_com: gen sh_sandy1 = sh_sandy if _n==1 // sandiness by judicial district
bys bailliage_id: egen bsandy = mean(sh_sandy1)
qui reghdfe nfert T01xaffected ruggedness [aw=wgt_4], nocons absorb(byear ${dist1} bailliage_id segment_affected_50) cluster(insee_com) nofoot
summ bsandy if e(sample)==1, d
local med = `r(p50)'
replace TxFxlarge = T01xaffected*(bsandy>`med')
replace TxFxsmall = T01xaffected*(bsandy<=`med')
* (1)
xi: reghdfe nfert $treatvars ruggedness sh_sandy  [aw=wgt_4], nocons absorb(byear ${dist1} bailliage_id segment_affected_50) cluster(insee_com) nofoot
test _b[TxFxlarge] = _b[TxFxsmall]
local pval = round(`r(p)',0.001)
summ nfert if e(sample)==1
local mdv = round(`r(mean)',0.01)
outreg2 using "$APPTAB/tableA9.tex", ///
    keep($treatvars ) ctitle(nfert) nocons tex(fragment) nonote nor2 noobs long dec(3) label replace ///
	addtext(p value difference, `pval', Cohort FE, Y, Border segment FE, Y, Judicial district FE, Y, Flexible trends, ., Soil conditions based on, Soil texture, Ruggedness, Y, Share sandy, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1, Mean dep var, `mdv') ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* (2)
xi: reghdfe nfert $treatvars ruggedness sh_sandy  [aw=wgt_5], nocons absorb(byear ${dist1} bailliage_id segment_affected_50 $controls) cluster(insee_com) nofoot
test _b[TxFxlarge] = _b[TxFxsmall]
local pval = round(`r(p)',0.001)
summ nfert if e(sample)==1
local mdv = round(`r(mean)',0.01)
outreg2 using "$APPTAB/tableA9.tex", ///
    keep($treatvars ) ctitle(nfert) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(p value difference, `pval', Cohort FE, Y, Border segment FE, Y, Judicial district FE, Y, Flexible trends, Y, Soil conditions based on, Soil texture, Ruggedness, Y, Share sandy, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1, Mean dep var, `mdv') ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
// based on ruggedness
qui reghdfe nfert T01xaffected ruggedness [aw=wgt_4], nocons absorb(byear ${dist1} bailliage_id segment_affected_50) cluster(insee_com) nofoot
summ ruggedness if e(sample)==1, d
local med = `r(p50)'
replace TxFxlarge = T01xaffected*(ruggedness<=`med')
replace TxFxsmall = T01xaffected*(ruggedness>`med')
* (3)
xi: reghdfe nfert $treatvars ruggedness sh_sandy [aw=wgt_4], nocons absorb(byear ${dist1} bailliage_id segment_affected_50) cluster(insee_com) nofoot
test _b[TxFxlarge] = _b[TxFxsmall]
local pval = round(`r(p)',0.001)
summ nfert if e(sample)==1
local mdv = round(`r(mean)',0.01)
outreg2 using "$APPTAB/tableA9.tex", ///
    keep($treatvars ) ctitle(nfert) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(p value difference, `pval', Cohort FE, Y, Border segment FE, Y, Judicial district FE, Y, Flexible trends, ., Soil conditions based on, Terrain ruggedness, Ruggedness, Y, Share sandy, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1, Mean dep var, `mdv') ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* (4)
xi: reghdfe nfert $treatvars ruggedness sh_sandy  [aw=wgt_5], nocons absorb(byear ${dist1} bailliage_id segment_affected_50 $controls) cluster(insee_com) nofoot
test _b[TxFxlarge] = _b[TxFxsmall]
local pval = round(`r(p)',0.001)
summ nfert if e(sample)==1
local mdv = round(`r(mean)',0.01)
outreg2 using "$APPTAB/tableA9.tex", ///
    keep($treatvars ) ctitle(nfert) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(p value difference, `pval', Cohort FE, Y, Border segment FE, Y, Judicial district FE, Y, Flexible trends, Y, Soil conditions based on, Terrain ruggedness, Ruggedness, Y, Share sandy, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1, Mean dep var, `mdv') ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ==============================================================================



********************************************************************************
*
* APPENDIX E. ROBUSTNESS OF RD-DD RESULTS, GENI DATABASE
*
********************************************************************************



* ==============================================================================
* Table E1: Covariate balance test for spatial regression-discontinuity analysis.
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* Table
label var affected "Reform"
local panels="A B"
local d = 17.49

forvalues o=1/2{

local pan `: word `o' of `panels''

* Bandwidths used in baseline RD regs
global da`o' = round(`d',0.01)
gen wgt_`o' = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_`o' = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_`o' = 0 if mi(wgt_`o')

* Balancedness RD regs
label var affected "RD estimate reform"
local t1 "wheat"
local s1 "price"
local t2 "refract"
local s2 "clergy"
local t3 "pre-reform"
local s3 "pdensity"
forvalues j=4/11{
local t`j' "distance"
}
local s4 "eveche"
local s5 "socpol"
local s6 "rebellion"
local s7 "bailliage"
local s8 "subdeleg"
local s9 "recette"
local s10 "post"
local s11 "road"
local re "replace"

* Wheat prices (individual level)
local j = 1
foreach y in logpwheat {
qui xi: reghdfe `y' affected [aw=wgt_`o'], absorb(byear ${dist`o'} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE1-panel`pan'.tex", ///
 	ctitle("`t`j''", "`s`j''") `re' ///
    keep(affected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Observations, `e(N)', N clusters, `e(N_clust)', Order polynomial, `o', Bandwidth, ${da`o'}, Cohort FE, Y, Border segment FE, Y, Unit of obs, Individual) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
local re " "
local j = `j'+1
}
* Municipality level variables
foreach y in cl_peril_i ldensity_1793 near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_post near_cassini {
preserve
collapse `y' affected wgt_`o' dista segment_affected_50 T01 dista_sq, by(insee_com)
label var affected "RD estimate reform"
qui xi: reghdfe `y' affected [aw=wgt_`o'], absorb(${dist`o'} segment_affected_50)
outreg2 using "$APPTAB/tableE1-panel`pan'.tex", ///
 	ctitle("`t`j''", "`s`j''") ///
    keep(affected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Observations, `e(N)', N clusters, ., Order polynomial, `o', Bandwidth, ${da`o'}, Cohort FE, ., Border segment FE, Y, Unit of obs, Locality) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) 
restore
local j = 1+`j'
}
local d = 30.53
}
* ==============================================================================




* ==============================================================================
* Table E2: Sensitivity to additional RD-DD specifications.
* ------------------------------------------------------------------------------

* ------------------------------------------------------------------------------
* Panel A. Two-dimensional running variable in latitude, longitude, latitude × longitude
* ------------------------------------------------------------------------------
use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* 2 dimensional polynomials
gen latlon = latitude*longitude
gen lat_sq = latitude^2
gen lon_sq = longitude^2
gen lat_sqlon = lat_sq*longitude
gen latlon_sq = latitude*lon_sq
global latlon1 c.latitude#T01#affected c.longitude#T01#affected c.latlon#T01#affected
global latlon2 c.latitude#T01#affected c.longitude#T01#affected c.lat_sq#T01#affected c.lon_sq#T01#affected c.latlon#T01#affected c.lat_sqlon#T01#affected c.latlon_sq#T01#affected  

* RD-DD regs
* (1)
qui reghdfe nfert affected T01xaffected if inrange(dista,-100,100), absorb(byear c.($latlon1 )#T01#affected) cluster(insee_com)
local nclust = `e(N_clust)'
local main = _b[affected]
margins , expression(_b[T01xaffected]/`main') post
outreg2 using "$APPTAB/tableE2-panelA.tex", ///
    tex(fragment) nonote nor2 noobs long dec(3) ///
	addtext(Cohort FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `nclust', Order polynomial, 1, Bandwidth, 100) ///
	addnote(Constant is ReformxFertile estimate; Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace
* (2)
qui reghdfe nfert affected T01xaffected if inrange(dista,-100,100), absorb(byear c.($latlon2 )#T01#affected) cluster(insee_com)
local nclust = `e(N_clust)'
local main2 = _b[affected]
margins, expression(_b[T01xaffected]/`main2') post
outreg2 using "$APPTAB/tableE2-panelA.tex", ///
    tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `nclust', Order polynomial, 2, Bandwidth, 100) ///
	addnote(Constant is ReformxFertile estimate; Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* (3)
qui reghdfe nfert affected T01xaffected if inrange(dista,-100,100), absorb(byear bailliage_id c.($latlon1 )#T01#affected) cluster(insee_com)
local nclust = `e(N_clust)'
margins, expression(_b[T01xaffected]/`main') post
outreg2 using "$APPTAB/tableE2-panelA.tex", ///
    tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `nclust', Order polynomial, 1, Bandwidth, 100) ///
	addnote(Constant is ReformxFertile estimate; Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* (4)
qui reghdfe nfert affected T01xaffected if inrange(dista,-100,100), absorb(byear $controls c.($latlon1 )#T01#affected) cluster(insee_com)
local nclust = `e(N_clust)'
local main4 = _b[affected]
margins, expression(_b[T01xaffected]/`main') post
qui summ nfert if e(sample)==1
local mdv = round(`r(mean)',0.01)
outreg2 using "$APPTAB/tableE2-panelA.tex", ///
    tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `nclust', Order polynomial, 1, Bandwidth, 100) ///
	addnote(Constant is ReformxFertile estimate; Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* (5)
qui reghdfe nfert affected T01xaffected if inrange(dista,-100,100), absorb(byear bailliage_id $controls c.($latlon1 )#T01#affected) cluster(insee_com)
local nclust = `e(N_clust)'
local main4 = _b[affected]
margins, expression(_b[T01xaffected]/`main') post
outreg2 using "$APPTAB/tableE2-panelA.tex", ///
    tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `nclust', Order polynomial, 1, Bandwidth, 100) ///
	addnote(Constant is ReformxFertile estimate; Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ------------------------------------------------------------------------------

	

* ------------------------------------------------------------------------------
* Panel B. Running variable in distance varies by years fertile post-reforms
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T#affected
global dist2 c.dista#T#affected c.dista_sq#T#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (2) 
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da2 = round(`d',0.01)
gen wgt_2 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
* col (3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da3 = round(`d',0.01)
gen wgt_3 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_3 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_3 = 0 if mi(wgt_3)

// baseline specifications + flexible trends
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
* col (4)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)
drop _A* _B* _C* _D*

* RD-DD results
* -------------------------------
* col (1) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_1], absorb(byear ${dist1} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelB.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace

* col (2) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_2], absorb(byear ${dist2} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelB.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 2, MSE Optimal bandwidth, $da2 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (3) 
xi: reghdfe nfert T01xaffected [aw=wgt_3], absorb(byear ${dist1} segment_affected_50 bailliage_id ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelB.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da3 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (4) 
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_4], absorb(byear ${dist1} segment_affected_50 $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelB.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da4 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
* col (5) 
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_5], absorb(byear ${dist1} segment_affected_50 bailliage_id $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelB.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da5 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ------------------------------------------------------------------------------


	
* ------------------------------------------------------------------------------
* Panel C. Uniform kernel
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(uniform) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1/`d' if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (2) 
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(uniform) bwselect(mserd)
local d = `e(h_mserd)'
global da2 = round(`d',0.01)
gen wgt_2 = 1/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1/`d' if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
* col (3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(uniform) bwselect(mserd)
local d = `e(h_mserd)'
global da3 = round(`d',0.01)
gen wgt_3 = 1/`d' if inrange(dista, 0, `d')
replace wgt_3 = 1/`d' if inrange(dista, -`d', 0)
replace wgt_3 = 0 if mi(wgt_3)

// baseline specifications + flexible trends
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
* col (4)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected logpwheat _A* _B* _C* _D*) kernel(uniform) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1/`d' if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(uniform) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1/`d' if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)
drop _A* _B* _C* _D*

* RD-DD results
* -------------------------------
label var T01xaffected "Reform X fertile post-reform"
* col (1) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_1], absorb(byear ${dist1} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelC.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace

* col (2) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_2], absorb(byear ${dist2} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelC.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 2, MSE Optimal bandwidth, $da2 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (3) 
xi: reghdfe nfert T01xaffected [aw=wgt_3], absorb(byear ${dist1} segment_affected_50 bailliage_id ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelC.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da3 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (4) 
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_4], absorb(byear ${dist1} segment_affected_50 $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelC.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da4 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
* col (5) 
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_5], absorb(byear ${dist1} segment_affected_50 bailliage_id $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelC.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da5 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ------------------------------------------------------------------------------


* ------------------------------------------------------------------------------
* Panel D. Henry sample - eighteenth-century cohorts
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
drop if (byear>1805) | (!inrange(byear,1705,1795) & birth_start_circa=="t") // for this robustness test
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (2) 
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da2 = round(`d',0.01)
gen wgt_2 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
* col (3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da3 = round(`d',0.01)
gen wgt_3 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_3 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_3 = 0 if mi(wgt_3)

// baseline specifications + flexible trends
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
* col (4)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(625)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(625)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)
drop _A* _B* _C* _D*


* RD-DD results
* -------------------------------
label var T01xaffected "Reform X fertile post-reform"
* col (1) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_1], absorb(byear ${dist1} segment_affected_50) cluster(insee_com)
qui summ nfert if e(sample)==1
outreg2 using "$APPTAB/tableE2-panelD.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace

* col (2) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_2], absorb(byear ${dist2} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelD.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 2, MSE Optimal bandwidth, $da2 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (3) 
xi: reghdfe nfert T01xaffected [aw=wgt_3], absorb(byear ${dist1} segment_affected_50 bailliage_id ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelD.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da3 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (4) 
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_4], absorb(byear ${dist1} segment_affected_50 $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelD.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da4 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
* col (5) 
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_5], absorb(byear ${dist1} segment_affected_50 bailliage_id $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panelD.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da5 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ------------------------------------------------------------------------------


	
* ------------------------------------------------------------------------------
* Panel E. Henry sample – rural municipalities (<20,000 inhabitants)
* Panel F. Henry sample – no arrondissement chief-lieux
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* Robustness cheks
global robustE = `"drop if pop_1793>=20000"'
global robustF = `"drop if rural_henry!=1"'

foreach pan in E F {
	
preserve

${robust`pan'}

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (2) 
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da2 = round(`d',0.01)
gen wgt_2 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
* col (3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da3 = round(`d',0.01)
gen wgt_3 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_3 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_3 = 0 if mi(wgt_3)

// baseline specifications + flexible trends
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
* col (4)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)
drop _A* _B* _C* _D*

* RD-DD results
* -------------------------------
label var T01xaffected "Reform X fertile post-reform"

* col (1) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_1], absorb(byear ${dist1} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panel`pan'.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace

* col (2) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_2], absorb(byear ${dist2} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panel`pan'.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 2, MSE Optimal bandwidth, $da2 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (3) 
xi: reghdfe nfert T01xaffected [aw=wgt_3], absorb(byear ${dist1} segment_affected_50 bailliage_id ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panel`pan'.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da3 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (4) 
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_4], absorb(byear ${dist1} segment_affected_50 $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panel`pan'.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da4 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
* col (5) 
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_5], absorb(byear ${dist1} segment_affected_50 bailliage_id $controls ) cluster(insee_com)
outreg2 using "$APPTAB/tableE2-panel`pan'.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da5 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

restore 
	
}
* ------------------------------------------------------------------------------

	
* ------------------------------------------------------------------------------
* Panel G. 100-kilometer border-segment fixed effects
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

replace segment_affected_50 = segment_affected_100

* sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
* -------------------------------
// baseline specifications
* col (1) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da1 = round(`d',0.01)
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
* col (2) 
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da2 = round(`d',0.01)
gen wgt_2 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
* col (3) 
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
global da3 = round(`d',0.01)
gen wgt_3 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_3 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_3 = 0 if mi(wgt_3)

// baseline specifications + flexible trends
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
* col (4)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* T01xaffected logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da4 = round(`d',0.01)
gen wgt_4 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_4 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_4 = 0 if mi(wgt_4)
* col (5)
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
global da5 = round(`d',0.01)
gen wgt_5 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_5 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_5 = 0 if mi(wgt_5)
drop _A* _B* _C* _D*

* RD-DD results
* -------------------------------
label var T01xaffected "Reform X fertile post-reform"

* col (1) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_1], absorb(byear ${dist1} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB\tableE2-panelG.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da1 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001) replace

* col (2) 
xi: reghdfe nfert affected T01xaffected [aw=wgt_2], absorb(byear ${dist2} segment_affected_50) cluster(insee_com)
outreg2 using "$APPTAB\tableE2-panelG.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 2, MSE Optimal bandwidth, $da2 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (3) 
xi: reghdfe nfert T01xaffected [aw=wgt_3], absorb(byear ${dist1} segment_affected_50 bailliage_id ) cluster(insee_com)
outreg2 using "$APPTAB\tableE2-panelG.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, ., N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da3 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)

* col (4) 
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_4], absorb(byear ${dist1} segment_affected_50 $controls ) cluster(insee_com)
outreg2 using "$APPTAB\tableE2-panelG.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, ., Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da4 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
	
* col (5) 
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_5], absorb(byear ${dist1} segment_affected_50 bailliage_id $controls ) cluster(insee_com)
outreg2 using "$APPTAB\tableE2-panelG.tex", ///
    keep(T01xaffected) nocons tex(fragment) nonote nor2 noobs long dec(3) label ///
	addtext(Cohort FE, Y, Border segment FE, Y, Bailliage FE, Y, Flexible trends, Y, N, `e(N)', N clusters, `e(N_clust)', Order polynomial, 1, MSE Optimal bandwidth, $da5 ) ///
	addnote(Geni h-sample women born 1700-1810 + geolocated, Triangular kernel; SE clustered by locality; *p<.05; **p<.01; ***p<.001)
* ==============================================================================



* ==============================================================================
* Figure E1: Balance RD plots.
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE bandwidth
global bw 17.49

* RD plots
keep if dista>-$bw & dista<$bw
gen disti = -dista
local t2 "Refractory clergy"
local d02 = 0
local d12 = 20
local d22 = 100
local ytit2 = "%"
local t3 "Pre-reform population density"
local d03 = 0
local d13 = 1
local d23 = 6
local ytit3 = "Log-population"
local t4 "Religious center within 15km"
local t5 "Political society within 15km"
local t6 "Rebellion 1779-1789 within 15km"
local t7 "Legal centre within 15km"
local t8 "Administrative center within 15km"
local t9 "Tax center within 15km"
local t10 "Paved road within 7.5km"
local t11 "Horse post within 15km"
forvalues t=4/11{
local d0`t' = 0
local d1`t' = 0.2
local d2`t' = 1
local ytit`t' = "share"
}

* RD BALANCE PLOTS
rdplot logpwheat disti, p(1) h($bw ) covs(segment_affected_50 baiFE_*) binselect(es) ///
	graph_options(xlabel(-15(5)15) ylabel(0.6(0.2)1.4) title("Local wheat price", color(black) size(huge)) ytitle("Log-price",  size(huge)) xtitle("Distance to inheritance border (> 0 reformed)", size(huge)) legend(off) graphr(lc(white) fc(white)) plotr(lc(black)) aspect(0.4))
graph export "$APPFIG/figureE1a.pdf", as(pdf) name("Graph") replace

local panels="a b c d e f g h i j k"
local j = 2
foreach y in cl_peril_i ldensity_1793 near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post {
local pan `: word `j' of `panels''
preserve 
collapse `y' disti segment_affected_50 baiFE_*, by(insee_com)
rdplot `y' disti, p(1) h($bw ) covs(segment_affected_50 baiFE_*) binselect(es) ///
	graph_options(xlabel(-15(5)15) ylabel(`d0`j''(`d1`j'')`d2`j'') title("`t`j''", color(black) size(huge)) ytitle("`ytit`j''", size(huge)) xtitle("Distance to inheritance border (> 0 reformed)", size(huge)) legend(off) graphr(lc(white) fc(white)) plotr(lc(black)) aspect(0.4))
graph export "$APPFIG/figureE1`pan'.pdf", as(pdf) name("Graph") replace
restore
local j = 1+`j'
}
* ==============================================================================



* ==============================================================================
* Figure E2: Trends in completed fertility under reformed and not reformed inheritance
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Base sample and vars for figure
keep if gender=="f" 
keep if byear>=1700 & byear<=1820
keep if twoplus_flag==1
gen T2 = T if impart==1
gen decade = 10 * floor(byear/10)

* Trends (all)
preserve
keep if dista!=.
bys decade: egen i1fert = mean(nfert) if affected==1
bys decade: egen i0fert = mean(nfert) if affected==0
collapse i1fert i0fert T2, by(decade byear)
keep if byear<=1810
twoway (scatter i0fert decade if byear>1704 & byear<=1810, m(o) mc(stc1)) ///
	   (scatter i1fert decade if byear>1704 & byear<=1810, m(o) mc(stc2)) ///
	   (lfit i0fert decade if decade<1760 & byear>1704, lc(stc1)) ///
	   (lfit i1fert decade if decade<1760 & byear>1704, lc(stc2)) ///
	   (line T2 byear, yaxis(2) lc(gray%25) lw(vthick)) ///
	   (pcarrowi 4.6 1754 4.6 1803 (1) " Cohorts fertile after reform", lcolor(black) mc(black) mlabcolor(black)), ///
	   ylabel(2.5(0.5)4.5, labs(medsmall) angle(horizontal)) ///
	   ytitle("Completed fertility", size(medlarge)) 	 ///
	   ytitle("Years fertile post-reform", axis(2) size(medlarge)) ylabel(0 "0" 5 "5" 10 "10" 15 "15" 20 "20" 25 "25", labs(small) angle(horizontal) axis(2)) ///
	   xlabel(1700(10)1810, labs(medsmall) angle(45)) ///
	   ymlabel(26 " ", axis(2) notick) ymlabel(4.6 " ", axis(1) notick) ///
	   title("(a) All", margin(large)) ///
	   legend(on order(2 "Reformed" 1 "Not reformed" 6 "Treatment (right)") cols(3) size(medsmall) pos(6) ring(1) region(lc(black))) ///
	   xline(1753.5, lcolor(black) lpattern(dash)) xtitle("Women's birth cohort", size(medium)) ///
	   aspectratio(0.5) graphregion(fcolor(white)) name(grall, replace)
restore	   
	   
* Trends (border)
preserve
keep if dista!=. & dista>-150 & dista<150
bys decade: egen i1fert = mean(nfert) if affected==1
bys decade: egen i0fert = mean(nfert) if affected==0
collapse i1fert i0fert T2, by(decade byear)
keep if byear<=1810
twoway (scatter i0fert decade if byear>1704 & byear<=1810, m(o) mc(stc1)) ///
	   (scatter i1fert decade if byear>1704 & byear<=1810, m(o) mc(stc2)) ///
	   (lfit i0fert decade if decade<1760 & byear>1704, lc(stc1)) ///
	   (lfit i1fert decade if decade<1760 & byear>1704, lc(stc2)) ///
	   (line T2 byear, yaxis(2) lc(gray%25) lw(vthick)) ///
	   (pcarrowi 4.6 1754 4.6 1803 (1) " Cohorts fertile after reform", lcolor(black) mc(black) mlabcolor(black)), ///
	   ylabel(2.5(0.5)4.5, labs(medsmall) angle(horizontal)) ///
	   ytitle("Completed fertility", size(medlarge)) 	 ///
	   ytitle("Years fertile post-reform", axis(2) size(medlarge)) ylabel(0 "0" 5 "5" 10 "10" 15 "15" 20 "20" 25 "25", labs(small) angle(horizontal) axis(2)) ///
	   xlabel(1700(10)1810, labs(medsmall) angle(45)) ///
	   ymlabel(26 " ", axis(2) notick) ymlabel(4.6 " ", axis(1) notick) ///
	   title("(b) 150 km of inheritance border", margin(large)) ///
	   legend(on order(2 "Reformed" 1 "Not reformed" 6 "Treatment (right)") cols(3) size(medsmall) pos(6) ring(1) region(lc(black))) ///
	   xline(1753.5, lcolor(black) lpattern(dash)) xtitle("Women's birth cohort", size(medium)) ///
	   aspectratio(0.5) graphregion(fcolor(white)) name(gr150, replace)

restore
grc1leg2 grall gr150, legendfrom(grall) name(grcomb, replace)
graph drop grall gr150
graph export "$APPFIG/figureE2.pdf", as(pdf) name("grcomb") replace
* ==============================================================================


* ==============================================================================
* Figure E3: Sensitivity to bandwidth choice for spatial RD-DD.
* ------------------------------------------------------------------------------

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

gen bw = . 

* (a) LINEAR POLYNOMIAL
* ------------------------------------------------------------------------------
local i = 1
forvalues d=-5(1)5{
qui{ 
	
local d1 = 19.47 + `d'
local d2 = 16.91 + `d'
local d3 = 20.6 + `d'
local d4 = 17.49 + `d'

forvalues j=1/4{
gen wgt_`j' = 1 - dista/`d`j'' if inrange(dista, 0, `d`j'')
replace wgt_`j' = 1 - dista/(-`d`j'') if inrange(dista, -`d`j'', 0)
replace wgt_`j' = 0 if mi(wgt_`j')
}

replace bw = `d' if _n==`i'

* (1) DD
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_1], absorb(byear $dist1 segment_affected_50 ) cluster(insee_com)
estimates store est1_`i'

* (2) DD + bailliage FE
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_2], absorb(byear $dist1 segment_affected_50 bailliage_id) cluster(insee_com)
estimates store est2_`i'

* (3) DD + flex trends
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_3], absorb(byear $dist1 segment_affected_50 $controls ) cluster(insee_com)
estimates store est3_`i'

* (4) DD + bailliage FE + flex trends
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_4], absorb(byear $dist1 segment_affected_50 bailliage_id $controls ) cluster(insee_com)
estimates store est4_`i'

drop wgt_*
local i = `i'+1
}
di `d'
}

* (a) QUADRATIC POLYNOMIAL
* ------------------------------------------------------------------------------
local i = 5
forvalues d=-10(2)10{
qui{ 
	
local d5 = 36.28 + `d'
local d6 = 29.29 + `d'
local d7 = 31.4 + `d'
local d8 = 30.53 + `d'

forvalues j=5/8{
gen wgt_`j' = 1 - dista/`d`j'' if inrange(dista, 0, `d`j'')
replace wgt_`j' = 1 - dista/(-`d`j'') if inrange(dista, -`d`j'', 0)
replace wgt_`j' = 0 if mi(wgt_`j')
}

replace bw = `d' if _n==`i'

* (5) DD
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_5], absorb(byear $dist2 segment_affected_50 ) cluster(insee_com)
estimates store est5_`i'

* (6) DD + bailliage FE
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_6], absorb(byear $dist2 segment_affected_50 bailliage_id) cluster(insee_com)
estimates store est6_`i'

* (7) DD + flex trends
xi: reghdfe nfert affected T01xaffected logpwheat [aw=wgt_7], absorb(byear $dist2 segment_affected_50 $controls ) cluster(insee_com)
estimates store est7_`i'

* (8) DD + bailliage FE + flex trends
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_8], absorb(byear $dist2 segment_affected_50 bailliage_id $controls ) cluster(insee_com)
estimates store est8_`i'

drop wgt_*
local i = `i'+1
}
di `d'
}


* graphs
local i = 1
forvalues d=-5(1)5{
global x`i' = 1+`d'*0.0835
local i = `i'+1
}

local title1 = "RD-DD"
local title2 = "RD-DD + bailliage FE"
local title3 = "RD-DD + flexible trends"
local title4 = "RD-DD + bailliage FE + flexible trends"
forvalues j=1/4{
coefplot (est`j'_1) (est`j'_2) (est`j'_3) (est`j'_4) (est`j'_5) ///
		 (est`j'_6, mcolor(stc2) msiz(medlarge) levels(90) ciopts(lc(stc2) lwidth(vvthin))) ///
		 (est`j'_7) (est`j'_8) (est`j'_9) (est`j'_10) (est`j'_11), ///
		 title("`title`j''") keep(T01xaffected) mcolor(stc1) msiz(medlarge) vertical ///
		 levels(90) ciopts(lc(stc1) lwidth(vthin)) /// 
		 xlabel($x1 "-5" $x2 "-4" $x3 "-3" $x4 "-2" $x5 "-1" $x6 "0" $x7 "+1" $x8 "+2" $x9 "+3" $x10 "+4" $x11 "+5") ylabel(-2.4(0.4)0.4) yline(0) ///
		 legend(on order(12 "MSE optimal bandwidth" 2 "-/+ km") rows(1) pos(6) ring(1) region(lc(black))) ///
		 name(g`j', replace) 
}
forvalues j=1/4{
local j2 = `j'+4
coefplot (est`j2'_5) (est`j2'_6) (est`j2'_7) (est`j2'_8) (est`j2'_9) ///
		 (est`j2'_10, mcolor(stc2) msiz(medlarge) levels(90) ciopts(lc(stc2) lwidth(vvthin))) ///
		 (est`j2'_11) (est`j2'_12) (est`j2'_13) (est`j2'_14) (est`j2'_15), ///
		 title("`title`j''") keep(T01xaffected) mcolor(stc1) msiz(medlarge) vertical ///
		 levels(90) ciopts(lc(stc1) lwidth(vthin)) /// 
		 xlabel($x1 "-10" $x2 "-8" $x3 "-6" $x4 "-4" $x5 "-2" $x6 "0" $x7 "+2" $x8 "+4" $x9 "+6" $x10 "+8" $x11 "+10") ylabel(-2.4(0.4)0.4) yline(0) ///
		 legend(on order(12 "MSE optimal bandwidth" 2 "-/+ km") rows(1) pos(6) ring(1) region(lc(black))) ///
		 name(g`j2', replace) 
}	

grc1leg2  g1 g2 g3 g4, rows(2) legendfrom(g1) name(a, replace)
graph export "$APPFIG/figureE3-a.pdf", as(pdf) name("a") replace

grc1leg2  g5 g6 g7 g8, rows(2) legendfrom(g5) name(b, replace)
graph export "$APPFIG/figureE3-b.pdf", as(pdf) name("b") replace
graph drop _all
* ==============================================================================


* ==============================================================================
* Figure E4: Conley adjusted standard errors with different distance cutoffs.
* ------------------------------------------------------------------------------

* ------------------------------------------------------------------------------
* (a) RD-DD linear
* ------------------------------------------------------------------------------

clear all
set maxvar 32767

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected
gen dpoly1 = dista if T01==0 & affected==0
gen dpoly2 = dista if T01==1 & affected==0
gen dpoly3 = dista if T01==0 & affected==1
gen dpoly4 = dista if T01==1 & affected==1
gen dpoly5 = dista^2 if T01==0 & affected==0
gen dpoly6 = dista^2 if T01==1 & affected==0
gen dpoly7 = dista^2 if T01==0 & affected==1
gen dpoly8 = dista^2 if T01==1 & affected==1
forvalues i=1/8{
	replace dpoly`i'=0 if dpoly`i'==.
}

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
forvalues o=1/2{
rdbwselect nfert dista, c(0) p(`o') vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_*) kernel(triangular) bwselect(mserd)
local d = `e(h_mserd)'
gen wgt_`o' = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_`o' = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_`o' = 0 if mi(wgt_`o')
}

* S.E. clustered by municipality
gen distcutoff = 0 if _n==1
// linear polynomial
xi: reghdfe nfert T01xaffected [aw=wgt_1], absorb(byear bailliage_id c.dista#T01#affected segment_affected_50) cluster(insee_com)
gen conleyz_linear = _b[T01xaffected]/_se[T01xaffected] if _n==1
// quadratic polynomial
xi: reghdfe nfert T01xaffected [aw=wgt_2], absorb(byear bailliage_id c.dista_sq#T01#affected c.dista#T01#affected segment_affected_50) cluster(insee_com)
gen conleyz_qua = _b[T01xaffected]/_se[T01xaffected] if _n==1

* Conley-adjusted SE, different cutoffs (over which spatial correlation assumed 0)
local n = 2
forvalues d=50(50)200{

di " ------------------------------ "
di " --------- cutoff `d' --------- "
di " ------------------------------ "

replace distcutoff = `d' if _n==`n'
// linear polynomial
acreg nfert T01xaffected dpoly1-dpoly4 baiFE_1-baiFE_136 baiFE_138-baiFE_379 [pw=wgt_1], spatial latitude(latitude) longitude(longitude) pfe1(byear) pfe2(segment_affected_50) id(profile_id) time(byear) lagcut(200) dist(`d')
replace conleyz_linear = _b[T01xaffected]/_se[T01xaffected] if _n==`n'
// quadratic polynomial
acreg nfert affected T01xaffected dpoly1-dpoly8 baiFE_1-baiFE_136 baiFE_138-baiFE_379 [pw=wgt_2], spatial latitude(latitude) longitude(longitude) pfe1(byear) pfe2(segment_affected_50) id(profile_id) time(byear) lagcut(200) dist(`d')
replace conleyz_qua = _b[T01xaffected]/_se[T01xaffected] if _n==`n'
local n = `n'+1
}

* save data 
keep distcutoff conleyz_linear conleyz_qua
keep if distcutoff!=.
save "$TEMP/RdDD1_conley-baiFE.dta", replace


* ------------------------------------------------------------------------------
* (b) RD-DD linear + flexible trends
* ------------------------------------------------------------------------------

clear all
set maxvar 32767

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected
gen dpoly1 = dista if T01==0 & affected==0
gen dpoly2 = dista if T01==1 & affected==0
gen dpoly3 = dista if T01==0 & affected==1
gen dpoly4 = dista if T01==1 & affected==1
gen dpoly5 = dista^2 if T01==0 & affected==0
gen dpoly6 = dista^2 if T01==1 & affected==0
gen dpoly7 = dista^2 if T01==0 & affected==1
gen dpoly8 = dista^2 if T01==1 & affected==1
forvalues i=1/8{
	replace dpoly`i'=0 if dpoly`i'==.
}

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)

* MSE optimal bandwidth & kernels
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
rdbwselect nfert dista, c(0) p(1) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
gen wgt_1 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_1 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_1 = 0 if mi(wgt_1)
drop _A* _B* _C* _D*

* S.E. clustered by municipality
gen distcutoff = 0 if _n==1
// linear polynomial
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_1], absorb(byear bailliage_id c.dista#T01#affected segment_affected_50 $controls ) cluster(insee_com)
gen conleyz_linear = _b[T01xaffected]/_se[T01xaffected] if _n==1
* save data 
preserve 
keep distcutoff conleyz_linear
keep if distcutoff!=.
save "$TEMP/RdDDflex_conley-baiFE.dta", replace
restore 

* Partial-out higher-dimension FE & flexible trends
encode bailliage_id, gen(nbailliage_id)
hdfe nfert T01xaffected logpwheat [aw=wgt_1], clear absorb(byear nbailliage_id c.dista#T01#affected segment_affected_50 $controls ) keepvars(byear profile_id insee_com wgt_1 distcutoff latitude longitude T01)
gen conleyz_linear = .

* Conley-adjusted SE, different cutoffs (over which spatial correlation assumed 0)
local n = 2
forvalues d=50(50)200{

di " ------------------------------ "
di " --------- cutoff `d' --------- "
di " ------------------------------ "

replace distcutoff = `d' if _n==`n'
// linear polynomial
acreg nfert T01xaffected logpwheat [pw=wgt_1], spatial latitude(latitude) longitude(longitude) id(profile_id) time(byear) lagcut(200) dist(`d')
replace conleyz_linear = _b[T01xaffected]/_se[T01xaffected] if _n==`n'
local n = `n'+1
}

* save data 
keep distcutoff conleyz_linear
keep if distcutoff!=.
append using "$TEMP/RdDDflex_conley-baiFE.dta"
save "$TEMP/RdDDflex_conley-baiFE.dta", replace


* ------------------------------------------------------------------------------
* (c) RD-DD quadratic
* ------------------------------------------------------------------------------

clear all
set maxvar 32767

use "$DAT/final-geni.dta", clear

* Sample
keep if gender=="f" 
keep if byear>=1700 & byear<=1810
keep if twoplus_flag==1

* geographic polynomials
global dist1 c.dista#T01#affected
global dist2 c.dista#T01#affected c.dista_sq#T01#affected
gen dpoly1 = dista if T01==0 & affected==0
gen dpoly2 = dista if T01==1 & affected==0
gen dpoly3 = dista if T01==0 & affected==1
gen dpoly4 = dista if T01==1 & affected==1
gen dpoly5 = dista^2 if T01==0 & affected==0
gen dpoly6 = dista^2 if T01==1 & affected==0
gen dpoly7 = dista^2 if T01==0 & affected==1
gen dpoly8 = dista^2 if T01==1 & affected==1
forvalues i=1/8{
	replace dpoly`i'=0 if dpoly`i'==.
}

* flexible trends
global controls c.(cl_peril_i pw1780 ldensity_1793)#byear near_eveche#byear near_socpol#byear near_rebellion#byear near_bailliage#byear near_subdeleg#byear near_recette#byear near_cassini#byear near_post#byear

* FE for bw calculation
qui encode insee_com, gen(ninsee_com)
qui tab byear, gen(byear_)
qui tab segment_affected_50, gen(segment_affected_50_)
qui tab bailliage_id, gen(baiFE_)


* MSE optimal bandwidth & kernels
xi i.byear*cl_peril_i, prefix(_A)
xi i.byear*ldensity_1793, prefix(_B)
xi i.byear*pw1780, prefix(_C)
local t= 1
foreach x in near_eveche near_socpol near_rebellion near_bailliage near_subdeleg near_recette near_cassini near_post{
xi i.`x'*i.byear, prefix(_D`t')
local t=`t'+1
}
rdbwselect nfert dista, c(0) p(2) vce(cluster ninsee_com) covs(byear_* segment_affected_50_* baiFE_* logpwheat _A* _B* _C* _D*) kernel(triangular) bwselect(mserd) bwcheck(650)
local d = `e(h_mserd)'
gen wgt_2 = 1 - dista/`d' if inrange(dista, 0, `d')
replace wgt_2 = 1 - dista/(-`d') if inrange(dista, -`d', 0)
replace wgt_2 = 0 if mi(wgt_2)
drop _A* _B* _C* _D*


* S.E. clustered by municipality
gen distcutoff = 0 if _n==1
// quadratic polynomial
xi: reghdfe nfert T01xaffected logpwheat [aw=wgt_2], absorb(byear bailliage_id c.dista_sq#T01#affected c.dista#T01#affected segment_affected_50 $controls ) cluster(insee_com)
gen conleyz_qua = _b[T01xaffected]/_se[T01xaffected] if _n==1
* save data 
preserve 
keep distcutoff conleyz_qua
keep if distcutoff!=.
save "$TEMP/RdDDflex_conley2-baiFE.dta", replace
restore 


* ------------------------------------------------------------------------------
* (d) RD-DD quadratic + flex. trends
* ------------------------------------------------------------------------------

* Partial-out higher-dimension FE & flexible trends
encode bailliage_id, gen(nbailliage_id)
hdfe nfert T01xaffected logpwheat [aw=wgt_2], clear absorb(byear nbailliage_id c.dista_sq#T01#affected c.dista#T01#affected segment_affected_50 $controls ) keepvars(byear profile_id insee_com wgt_2 distcutoff latitude longitude T01)
gen conleyz_qua = .

* Conley-adjusted SE, different cutoffs (over which spatial correlation assumed 0)
local n = 2
forvalues d=50(50)200{

di " ------------------------------ "
di " --------- cutoff `d' --------- "
di " ------------------------------ "

replace distcutoff = `d' if _n==`n'
// linear polynomial
acreg nfert T01xaffected logpwheat [pw=wgt_2], spatial latitude(latitude) longitude(longitude) id(profile_id) time(byear) lagcut(200) dist(`d')
replace conleyz_qua = _b[T01xaffected]/_se[T01xaffected] if _n==`n'
local n = `n'+1
}
* save data 
preserve 
keep distcutoff conleyz_qua
keep if distcutoff!=.
append using "$TEMP/RdDDflex_conley2-baiFE.dta"
save "$TEMP/RdDDflex_conley2-baiFE.dta", replace
restore 


* ------------------------------------------------------------------------------
* Graphs
* ------------------------------------------------------------------------------

// linear
use "$TEMP/RdDD1_conley-baiFE.dta", clear
twoway (scatter conleyz_linear distcutoff, msiz(medium)) (scatter conleyz_linear distcutoff if distcutoff==0, msiz(medium)), ///
	xtitle(Distance cutoff) ytitle(Z (coef/s.e.)) ylabel(-3(1)0) ///
	xlabel(0(50)200) xmlabel(-10 " " 210 " ", notick) scale(1.5) ///
	yline(-1.645, lcolor(gs10)) ///
	legend(on order(2 "Baseline SE (clustered by municipality)" 1 "Conley-adjusted SE") rows(2) ring(0) pos(2) region(lc(black)) siz(medsmall)) name(a, replace)
graph export "$APPFIG/figureE4a.png", replace

// linear + flex trends
use "$TEMP/RdDDflex_conley-baiFE.dta", clear
twoway (scatter conleyz_linear distcutoff, msiz(medium)) (scatter conleyz_linear distcutoff if distcutoff==0, msiz(medium)), ///
	xtitle(Distance cutoff) ytitle(Z (coef/s.e.)) ylabel(-3(1)0) ///
	xlabel(0(50)200) xmlabel(-10 " " 210 " ", notick) scale(1.5) ///
	yline(-1.645, lcolor(gs10)) ///
	legend(on order(2 "Baseline SE (clustered by municipality)" 1 "Conley-adjusted SE") rows(2) ring(0) pos(2) region(lc(black)) siz(medsmall)) name(b, replace)
graph export "$APPFIG/figureE4b.png", replace

// quadratic
use "$TEMP/RdDD1_conley-baiFE.dta", clear
twoway (scatter conleyz_qua distcutoff, msiz(medium)) (scatter conleyz_qua distcutoff if distcutoff==0, msiz(medium)), ///
	xtitle(Distance cutoff) ytitle(Z (coef/s.e.)) ylabel(-3(1)0) ///
	xlabel(0(50)200) xmlabel(-10 " " 210 " ", notick) scale(1.5) ///
	yline(-1.645, lcolor(gs10)) ///
	legend(on order(2 "Baseline SE (clustered by municipality)" 1 "Conley-adjusted SE") rows(2) ring(0) pos(2) region(lc(black)) siz(medsmall)) name(c, replace)
graph export "$APPFIG/figureE4c.png", replace

// quadratic + flex trends
use "$TEMP/RdDDflex_conley2-baiFE.dta", clear
twoway (scatter conleyz_qua distcutoff, msiz(medium)) (scatter conleyz_qua distcutoff if distcutoff==0, msiz(medium)), ///
	xtitle(Distance cutoff) ytitle(Z (coef/s.e.)) ylabel(-3(1)0) ///
	xlabel(0(50)200) xmlabel(-10 " " 210 " ", notick) scale(1.5) ///
	yline(-1.645, lcolor(gs10)) ///
	legend(on order(2 "Baseline SE (clustered by municipality)" 1 "Conley-adjusted SE") rows(2) ring(0) pos(2) region(lc(black)) siz(medsmall)) name(d, replace)
graph export "$APPFIG/figureE4d.png", replace

graph drop _all
* ==============================================================================

timer off 1
timer list